Analysing the effect of treatment on molecular scores using aligned-rank transformed analysis of variance

1 Analysing the effect of treatment on molecular scores using aligned-rank transformed analysis of variance

1.1 INTRODUCTION

This analysis assesses how felzartamab therapy affects molecular sores measured with gene expression data from biopsies from kidney allografts taken from patient with antibody mediated rejection. The assumptions of normality and heterogeneity of variance differ across molecular scores, thus we avail to alinged-rank transformation to carry out a non-parametric ANOVA. This is powerful in this context because it allows us to test for the interactive effect of treatment (i.e., placebo vs felzartamab) and time (i.e., change in scores from biopsy to biopsy).

References:

• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).

• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).

• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).

1.2 HOUSEKEEPING

• First we need to load the necessary R packages and data.

# Load necessary libraries
library(tidyverse)
library(broom)
library(rstatix)
library(flextable)
library(officer)
library(emmeans)
library(multcomp)
library(ARTool)
library(rcompanion)
library(Biobase)
library(ggpubr) 
library(patchwork) 
library(ggprism) 

# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0

# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)

# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")

# Seed for reproducibility
set.seed(42)

• Now we define the categories of molecular scores for easier organization downstream.

# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")

# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)

1.3 DATA WRANGLING

• Now we process the data so we can iteratively carry out the analyses on each score.

# DEFINE THE DATA ####
data <- set %>%
    pData()  %>% 
    dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
    dplyr::filter(Patient %nin% c(15, 18)) %>%
    left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
    relocate(sample_pairs, .after = Felzartamab_Followup) %>%
    arrange(Felzartamab, Patient, Group)


# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
    expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
    nest(.by = category) %>%
    mutate(
        features = map(
            category,
            function(category) {
                if (category == "cfDNA") {
                    vars_cfDNA
                } else if (category == "ABMR") {
                    vars_abmr
                } else if (category == "TCMR") {
                    vars_tcmr
                } else if (category == "macrophage") {
                    vars_macrophage
                } else if (category == "injury") {
                    vars_injury
                } else if (category == "parenchyma") {
                    vars_parenchyma
                } 
            }
        ),
        data = pmap(
            list(features, data),
            function(features, data) {
                data %>%
                    dplyr::select(
                        Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
                        all_of(features)
                    )
            }
        )
    )


# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
    mutate(
        data_univariate = pmap(
            list(data, features),
            function(data, features) {
                data %>%
                    pivot_longer(
                        cols = all_of(features),
                        names_to = "variable",
                        values_to = "value"
                    ) %>%
                    nest(.by = variable) %>%
                    mutate(
                        variable = variable %>%
                            factor(
                                levels = vars
                            ),
                        annotation = case_when(
                            variable %in% vars_cfDNA ~ "cfDNA",
                            variable %in% vars_tcmr ~ "TCMR-related",
                            variable %in% vars_abmr ~ "ABMR-related",
                            variable %in% vars_macrophage ~ "macrophage-related",
                            variable %in% vars_injury ~ "injury-related",
                            variable %in% vars_parenchyma ~ "parenchyma-related",
                            TRUE ~ " "
                        ) %>%
                            factor(
                                levels = c(
                                    "cfDNA",
                                    "ABMR-related",
                                    "TCMR-related",
                                    "macrophage-related",
                                    "injury-related",
                                    "parenchyma-related",
                                    "archetypes",
                                    "rejection PC",
                                    "injury PC"
                                )
                            ),
                        score = case_when(
                            variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
                            variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
                            variable == "TCB" ~ "T cell burden (TCB)",
                            variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
                            variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
                            variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
                            variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
                            variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
                            variable == "NKB" ~ "NK cell burden (NKB)",
                            variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
                            variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
                            variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
                            variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
                            variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
                            variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
                            variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
                            variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
                            variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
                            variable == "KT1" ~ "Kidney parenchymal (KT1)",
                            variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
                        ), .before = 1
                    ) %>%
                    arrange(annotation, variable)
            }
        )
    ) %>%
    dplyr::select(category, data_univariate) %>%
    unnest(data_univariate) %>%
    mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)

• This is what the data look like after we have organized them. Each row of the dataframe contains the data for each score, the category it belongs to, it’s category annotation, it’s long-format name, and the name of the variable in the data.

options(pillar.max_columns = NULL)
options(pillar.width = NULL)
data_01 %>% print(n="all")

## # A tibble: 20 x 6
##    category   annotation         score variable n_cat data    
##    <chr>      <fct>              <chr> <fct>    <int> <list>  
##  1 cfDNA      cfDNA              Dono~ cfDNA_c~     1 <tibble>
##  2 ABMR       ABMR-related       ABMR~ ABMRpm       5 <tibble>
##  3 ABMR       ABMR-related       Glom~ ggt0         5 <tibble>
##  4 ABMR       ABMR-related       Peri~ ptcgt0       5 <tibble>
##  5 ABMR       ABMR-related       DSA-~ DSAST        5 <tibble>
##  6 ABMR       ABMR-related       NK c~ NKB          5 <tibble>
##  7 TCMR       TCMR-related       TCMR~ TCMRt        5 <tibble>
##  8 TCMR       TCMR-related       Tubu~ tgt1         5 <tibble>
##  9 TCMR       TCMR-related       Inte~ igt1         5 <tibble>
## 10 TCMR       TCMR-related       Cyto~ QCAT         5 <tibble>
## 11 TCMR       TCMR-related       T ce~ TCB          5 <tibble>
## 12 macrophage macrophage-related Alte~ AMAT1        2 <tibble>
## 13 macrophage macrophage-related Cons~ QCMAT        2 <tibble>
## 14 injury     injury-related     Inju~ IRRAT30      5 <tibble>
## 15 injury     injury-related     Inju~ IRITD3       5 <tibble>
## 16 injury     injury-related     Inju~ IRITD5       5 <tibble>
## 17 injury     injury-related     Fibr~ cigt1        5 <tibble>
## 18 injury     injury-related     Atro~ ctgt1        5 <tibble>
## 19 parenchyma parenchyma-related Kidn~ KT1          2 <tibble>
## 20 parenchyma parenchyma-related Kidn~ KT2          2 <tibble>

1.4 CALCULATE MEDIAN SCORES AND MEDIAN EFFECT OF TREATMENT

• At this stage we can summarize the median values for each score by treatment and follow-up.

# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
    mutate(
        medians = map(
            data,
            function(data) {
                data %>%
                    reframe(
                        median = value %>% median(),
                        IQR = value %>% IQR(),
                        sample_pairs = n(),
                        .by = c(Followup, Felzartamab, Felzartamab_Followup)
                    ) %>%
                    arrange(Felzartamab_Followup)
            }
        ),
        medians_delta = map(
            medians,
            function(medians) {
                medians %>%
                    reframe(
                        median_delta = combn(median, 2, diff) %>% as.numeric(),
                        IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
                        sample_pairs = sample_pairs,
                        .by = c(Felzartamab)
                    ) %>%
                    mutate(
                        Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
                            factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
                        .before = sample_pairs
                    ) %>%
                    mutate(
                        median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
                        IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
                        .by = c(Followup_pairwise),
                        .before = sample_pairs
                    )
            }
        )
    )

• Let’s take a look at what the median summaries look like for the ABMRprob score.

data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians)

## [[1]]
## # A tibble: 6 x 6
##   Followup Felzartamab Felzartamab_Followup median   IQR
##   <fct>    <fct>       <fct>                 <dbl> <dbl>
## 1 Baseline Placebo     Baseline_Placebo      0.670 0.355
## 2 Week24   Placebo     Week24_Placebo        0.774 0.390
## 3 Week52   Placebo     Week52_Placebo        0.518 0.498
## 4 Baseline Felzartamab Baseline_Felzartamab  0.736 0.308
## 5 Week24   Felzartamab Week24_Felzartamab    0.167 0.357
## 6 Week52   Felzartamab Week52_Felzartamab    0.702 0.539
## # i 1 more variable: sample_pairs <int>

• Now we can look at the change in median ABMRprob scores in the placebo and felzartamab arms, as well as the differnce in those medians (i.e., ΔΔ median - this is our summarized treatment effect).

options(pillar.max_columns = Inf)
options(pillar.width = Inf)
data_02 %>%
    dplyr::filter(variable == "ABMRpm") %>%
    pull(medians_delta)

## [[1]]
## # A tibble: 6 x 7
##   Felzartamab median_delta IQR_delta Followup_pairwise
##   <fct>              <dbl>     <dbl> <fct>            
## 1 Placebo           0.105      0.373 Baseline - Week24
## 2 Placebo          -0.151      0.427 Baseline - Week52
## 3 Placebo          -0.256      0.444 Week24 - Week52  
## 4 Felzartamab      -0.568      0.332 Baseline - Week24
## 5 Felzartamab      -0.0335     0.423 Baseline - Week52
## 6 Felzartamab       0.535      0.448 Week24 - Week52  
##   median_delta_delta IQR_delta_delta sample_pairs
##                <dbl>           <dbl>        <int>
## 1             -0.673           0.353           10
## 2              0.118           0.425           10
## 3              0.791           0.446           10
## 4             -0.673           0.353           10
## 5              0.118           0.425           10
## 6              0.791           0.446           10

1.5 CARRY OUT ALIGNED-RANK TRANSFORM ANOVA

• Now we carry out the analysis. This includes testing the specific contrasts for the interaction of treatment and time (i.e., our effect of treatment).

# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
    mutate(
        art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
        art_aov = map(art, anova, type = "II"),
        art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
        art_aov_contrast = map(
            art,
            art.con,
            formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
        ),
        art_aov_contrast_tidy = map(art_aov_contrast, tidy),
        art_aov_contrast_cld = map(
            art_aov_contrast_tidy,
            function(art_aov_contrast_tidy) {
                art_aov_contrast_tidy %>%
                    as.data.frame() %>%
                    cldList(adj.p.value ~ Followup:Felzartamab, data = .)
            }
        )
    )

1.6 CREATE A SUMMARY TABLE OF ANOVA RESULTS FOR EACH SCORE

# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
    dplyr::select(annotation, n_cat, score, art_aov) %>%
    unnest(everything()) %>%
    dplyr::rename(p.value = `Pr(>F)`) %>%
    mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
    dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
    flextable::flextable() %>%
    flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
    flextable::merge_h(part = "header") %>%
    flextable::merge_v(j = 1:2) %>%
    flextable::fontsize(size = 8, part = "all") %>%
    flextable::align(align = "center", part = "all") %>%
    flextable::bg(bg = "white", part = "all") %>%
    flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
    flextable::colformat_double(j = 2:3, digits = 2) %>%
    flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
    flextable::border_remove() %>%
    flextable::bold(part = "header") %>%
    flextable::padding(padding = 0, part = "all") %>%
    flextable::border(border = officer::fp_border(), part = "all")  %>%
    flextable::autofit()

Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients
annotation	score	Term	F	p.value	FDR
cfDNA	Donor-derived cell-free DNA (dd-cfDNA, cp/mL)	Followup	4.418	0.019	0.029
		Felzartamab	1.598	0.222	0.222
		Followup:Felzartamab	10.396	0.000	0.001
ABMR-related	ABMR classifier (ABMRProb)	Followup	6.659	0.003	0.007
		Felzartamab	1.620	0.219	0.235
		Followup:Felzartamab	9.605	0.000	0.001
	Glomerulitis classifier (g>0Prob)	Followup	11.768	0.000	0.001
		Felzartamab	1.511	0.235	0.235
		Followup:Felzartamab	10.743	0.000	0.001
	Peritubular capillaritis classifier (ptc>0Prob)	Followup	14.188	0.000	0.000
		Felzartamab	3.178	0.092	0.125
		Followup:Felzartamab	7.749	0.002	0.004
	DSA-selective transcripts (DSAST)	Followup	10.508	0.000	0.001
		Felzartamab	1.747	0.203	0.234
		Followup:Felzartamab	4.179	0.023	0.035
	NK cell burden (NKB)	Followup	6.512	0.004	0.007
		Felzartamab	1.956	0.179	0.224
		Followup:Felzartamab	5.114	0.011	0.018
TCMR-related	TCMR classifier (TCMRProb)	Followup	3.555	0.039	0.165
		Felzartamab	1.417	0.249	0.468
		Followup:Felzartamab	1.124	0.336	0.504
	Tubulitis classifier (t>1Prob)	Followup	1.016	0.372	0.506
		Felzartamab	1.853	0.190	0.408
		Followup:Felzartamab	0.065	0.937	0.968
	Interstitial infiltrate classifier (i>1Prob)	Followup	3.414	0.044	0.165
		Felzartamab	0.997	0.331	0.504
		Followup:Felzartamab	0.032	0.968	0.968
	Cytotoxic T cell-associated (QCAT)	Followup	5.693	0.007	0.107
		Felzartamab	4.118	0.057	0.172
		Followup:Felzartamab	0.927	0.405	0.506
	T cell burden (TCB)	Followup	4.360	0.020	0.151
		Felzartamab	2.868	0.108	0.269
		Followup:Felzartamab	0.301	0.742	0.856
macrophage-related	Alternatively activated macrophage (AMAT1)	Followup	1.549	0.226	0.304
		Felzartamab	1.047	0.320	0.320
		Followup:Felzartamab	2.406	0.105	0.209
	Constitutive macrophage (QCMAT)	Followup	3.378	0.045	0.209
		Felzartamab	3.100	0.095	0.209
		Followup:Felzartamab	1.426	0.253	0.304
injury-related	Injury-repair associated (IRRAT30)	Followup	0.055	0.946	0.946
		Felzartamab	3.842	0.066	0.223
		Followup:Felzartamab	1.915	0.162	0.347
	Injury-repair induced, day 3 (IRITD3)	Followup	0.174	0.841	0.946
		Felzartamab	3.268	0.087	0.223
		Followup:Felzartamab	3.125	0.056	0.223
	Injury-repair induced, day 5 (IRITD5)	Followup	0.073	0.930	0.946
		Felzartamab	1.088	0.311	0.583
		Followup:Felzartamab	2.586	0.089	0.223
	Fibrosis classifier (ci>1Prob)	Followup	0.360	0.700	0.876
		Felzartamab	3.482	0.078	0.223
		Followup:Felzartamab	0.938	0.401	0.668
	Atrophy classifier (ct>1Prob)	Followup	0.568	0.572	0.780
		Felzartamab	3.297	0.086	0.223
		Followup:Felzartamab	0.602	0.553	0.780
parenchyma-related	Kidney parenchymal (KT1)	Followup	0.968	0.389	0.518
		Felzartamab	1.628	0.218	0.518
		Followup:Felzartamab	0.605	0.552	0.552
	Kidney parenchymal - no solute carriers (KT2)	Followup	0.859	0.432	0.518
		Felzartamab	0.911	0.352	0.518
		Followup:Felzartamab	1.038	0.365	0.518

1.7 SUMMARISE THE EFFECT OF TREATMENT FOR EACH SCORE ACROSS ALL TIME PERIODS

• Here we present the effect of treatment (ΔΔ) for each score, along with the effect of time (Δ) for placebo and felzartmab arms.

flextable_pairwise

Table 1. Effect of felzartamab on molecular scores
Annotation	Score	Baseline - Week 24				Week 24 - Week 52				Baseline - Week 52
		Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR	Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR	Δ Placebo (N=10)	Δ Felzartamab (N=10)	ΔΔ	ΔΔ FDR
cfDNA	Donor-derived cell-free DNA (dd-cfDNA, cp/mL)	0.50 (76.12)	-55.00 (30)	-55.50 (53.06)	2e-04	-3.50 (69.12)	14.50 (25.5)	18.00 (47.31)	0.027	-3.00 (51)	-40.50 (45.75)	-37.50 (48.38)	0.045
ABMR-related	ABMR classifier (ABMRProb)	0.10 (0.37)	-0.57 (0.33)	-0.67 (0.35)	1e-03	-0.26 (0.44)	0.54 (0.45)	0.79 (0.45)	1e-03	-0.15 (0.43)	-0.03 (0.42)	0.12 (0.43)	0.709
	Glomerulitis classifier (g>0Prob)	-0.02 (0.23)	-0.31 (0.35)	-0.29 (0.29)	3e-03	-0.19 (0.34)	0.31 (0.41)	0.50 (0.37)	2e-04	-0.22 (0.3)	0.00 (0.38)	0.21 (0.34)	0.277
	Peritubular capillaritis classifier (ptc>0Prob)	-0.12 (0.25)	-0.55 (0.27)	-0.43 (0.26)	3e-03	-0.10 (0.37)	0.37 (0.31)	0.47 (0.34)	3e-03	-0.22 (0.3)	-0.18 (0.35)	0.04 (0.33)	0.977
	DSA-selective transcripts (DSAST)	-0.09 (0.22)	-0.37 (0.24)	-0.28 (0.23)	0.026	-0.16 (0.32)	0.16 (0.33)	0.32 (0.33)	0.026	-0.24 (0.33)	-0.20 (0.25)	0.04 (0.29)	0.988
	NK cell burden (NKB)	0.02 (0.22)	-0.78 (0.49)	-0.80 (0.36)	0.014	-0.19 (0.4)	0.68 (0.59)	0.87 (0.49)	0.014	-0.17 (0.46)	-0.10 (0.62)	0.07 (0.54)	0.972
TCMR-related	TCMR classifier (TCMRProb)	-0.01 (0.02)	0.00 (0)	0.01 (0.01)	0.358	0.00 (0.02)	0.00 (0)	0.00 (0.01)	0.856	-0.01 (0.02)	0.00 (0)	0.01 (0.01)	0.358
	Tubulitis classifier (t>1Prob)	-0.01 (0.06)	-0.02 (0.04)	-0.01 (0.05)	1.000	0.00 (0.04)	0.01 (0.02)	0.01 (0.03)	1.000	-0.01 (0.05)	-0.01 (0.03)	-0.01 (0.04)	1.000
	Interstitial infiltrate classifier (i>1Prob)	-0.01 (0.04)	-0.02 (0.03)	-0.01 (0.04)	0.964	0.01 (0.04)	0.01 (0.04)	0.00 (0.04)	0.964	0.00 (0.05)	-0.01 (0.06)	-0.01 (0.05)	0.964
	Cytotoxic T cell-associated (QCAT)	-0.20 (0.41)	-0.70 (0.46)	-0.50 (0.43)	0.425	0.10 (0.33)	0.54 (0.57)	0.44 (0.45)	0.425	-0.10 (0.46)	-0.16 (0.41)	-0.06 (0.44)	0.871
	T cell burden (TCB)	-0.44 (0.71)	-0.79 (0.56)	-0.36 (0.64)	0.825	0.28 (0.57)	0.34 (0.6)	0.06 (0.58)	0.905	-0.16 (0.88)	-0.45 (0.43)	-0.29 (0.65)	0.825
macrophage-related	Alternatively activated macrophage (AMAT1)	0.08 (0.2)	-0.27 (0.45)	-0.35 (0.33)	0.128	0.09 (0.23)	0.11 (0.47)	0.02 (0.35)	0.812	0.18 (0.19)	-0.16 (0.32)	-0.33 (0.25)	0.128
	Constitutive macrophage (QCMAT)	0.00 (0.15)	-0.24 (0.27)	-0.24 (0.21)	0.340	0.01 (0.12)	0.09 (0.27)	0.08 (0.2)	0.685	0.01 (0.11)	-0.15 (0.24)	-0.16 (0.18)	0.348
injury-related	Injury-repair associated (IRRAT30)	0.18 (0.52)	0.05 (0.4)	-0.13 (0.46)	0.353	0.05 (0.66)	-0.11 (0.72)	-0.16 (0.69)	0.353	0.23 (0.6)	-0.06 (0.48)	-0.29 (0.54)	0.175
	Injury-repair induced, day 3 (IRITD3)	0.03 (0.19)	0.04 (0.12)	0.01 (0.15)	0.945	0.03 (0.18)	-0.18 (0.18)	-0.20 (0.18)	0.060	0.06 (0.19)	-0.13 (0.14)	-0.19 (0.16)	0.060
	Injury-repair induced, day 5 (IRITD5)	-0.01 (0.17)	0.04 (0.14)	0.05 (0.15)	0.932	0.07 (0.2)	-0.17 (0.21)	-0.24 (0.21)	0.093	0.06 (0.18)	-0.13 (0.17)	-0.18 (0.18)	0.093
	Fibrosis classifier (ci>1Prob)	-0.02 (0.29)	0.13 (0.32)	0.15 (0.3)	0.434	0.06 (0.25)	-0.08 (0.53)	-0.14 (0.39)	0.844	0.04 (0.25)	0.05 (0.37)	0.01 (0.31)	0.434
	Atrophy classifier (ct>1Prob)	0.01 (0.33)	0.12 (0.31)	0.12 (0.32)	0.587	0.09 (0.27)	-0.06 (0.44)	-0.15 (0.36)	0.587	0.10 (0.32)	0.07 (0.41)	-0.03 (0.36)	0.587
parenchyma-related	Kidney parenchymal (KT1)	-0.05 (0.22)	-0.15 (0.23)	-0.09 (0.22)	0.731	-0.01 (0.11)	0.10 (0.22)	0.11 (0.17)	0.705	-0.06 (0.25)	-0.04 (0.18)	0.02 (0.22)	0.705
	Kidney parenchymal - no solute carriers (KT2)	-0.09 (0.39)	-0.16 (0.37)	-0.07 (0.38)	0.577	-0.09 (0.2)	0.14 (0.38)	0.23 (0.29)	0.484	-0.18 (0.42)	-0.02 (0.31)	0.16 (0.37)	0.577
Grey shading denotes ANOVA interactive effect FDR < 0.05 FDR correction was carried out within each annotation grouping

1.8 PLOT THE RESULTS

# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
    dplyr::filter(variable == "AMAT1") %>%
    pull(plot_violin) %>%
    ggpubr::get_legend() %>%
    ggpubr::as_ggplot() +
    theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))


# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("ABMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
        plot.background = element_rect(fill = "grey95", colour = " white")
    )

panel_violin_abmr <- panel_violin_abmr %>%
    ggarrange(
        labels = "B",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
    dplyr::filter(category %in% c("TCMR")) %>%
    pull(plot_violin) %>%
    wrap_plots(nrow = 1, ncol = 5) &
    theme(
        axis.text = element_text(size = 10, colour = "black"),
        legend.position = "none",
    )

panel_violin_tcmr <- panel_violin_tcmr %>%
    ggarrange(
        labels = "C",
        font.label = list(size = 25, face = "bold"),
        legend = "none"
    ) %>%
    ggpubr::annotate_figure(
        top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
            face = "bold.italic",
            size = 25,
            hjust = 1.27
        )
    )


# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
    panel_violin_abmr,
    panel_violin_tcmr,
    nrow = 2,
    heights = c(1, 1)
)

ggarrange(
    panel_legend,
    panels_violin,
    nrow = 2,
    heights = c(0.125, 1)
)